import numpy as np
import matplotlib.pyplot as plt
# Returns the num channels for a passed in image np ndarray
def get_num_channels(img):
# Check for grayscale
if len(img.shape) < 3:
return 1
return img.shape[-1]
# Convert images to grayscale and print shape
def convert_grayscale(image):
weights = [0.3, 0.6, 0.1]
channels = get_num_channels(image)
if channels == 1: # Grayscale
return image
elif channels == 3: # RGB
return (image @ weights) / 255.0
elif channels == 4: # ARGB
# Only use the first 3 channels
return (image[:,:,:3] @ weights) /255.0
PATH = './cell_images'
files = [
'0001.000.png',
'0001.001.png',
'0001.002.png',
'0001.004.png',
'0001.005.png',
'0001.006.png',
]
images = []
for fname in files:
image = plt.imread(f'{PATH}/{fname}')
print(f'{fname}')
print(f'shape: {image.shape}')
print(f'max: {image.max()}')
print(f'min: {image.min()}')
print()
# Convert images to grayscale
images.append(convert_grayscale(image))
0001.000.png shape: (510, 510) max: 1.0 min: 0.0 0001.001.png shape: (510, 510) max: 1.0 min: 0.0 0001.002.png shape: (510, 510) max: 1.0 min: 0.0 0001.004.png shape: (510, 510) max: 1.0 min: 0.0 0001.005.png shape: (510, 510) max: 1.0 min: 0.0 0001.006.png shape: (510, 510) max: 1.0 min: 0.0
def display_images(imgs, fig_name):
rows = int(np.ceil(len(images) / 4))
fig, axs = plt.subplots(rows,4, figsize=(rows*4,4))
fig.tight_layout()
for i, axe in enumerate(axs):
for j, ax in enumerate(axe):
idx = (i*4) + j
if idx < len(imgs):
ax.imshow(imgs[idx], cmap='gray')
ax.set_title(files[idx])
# ax.set_xlabel(files[idx])
plt.suptitle('Cell Images')
plt.tight_layout()
plt.savefig(fig_name)
display_images(images, 'gs_images.png')
Read in and create overlapping crops and images from the drone image as another example.
def create_crops_from_image(img_path):
oth_images = []
oimage = plt.imread(img_path)
for i in range(10):
x = np.random.randint(2000, 2400)
y = np.random.randint(1000, 1400)
oth_images.append(oimage[x:x+500, y:y+500])
plt.imsave(f'./other_images/drone_00{i}.png', oth_images[i])
# Convert images to grayscale to use
gs_images = []
for i in range(len(oth_images)):
gs_images.append(convert_grayscale(oth_images[i]))
if i < 3:
plt.figure()
plt.imshow(gs_images[i], cmap='gray')
return gs_images
other_images = create_crops_from_image('./other_images/drone.jpg')
I followed the book for how to do the low pass filter in the Fourier domain and display multiple images power spectrums, the filters, and the resulting image back in the spatial domain after it has been filtered in the Fourier domain. I implemented the main 3 types of low pass filters described in the textbook which are the ideal low pass filter (ILPF), gaussian low pass filter (GLPF), and Butterworth low pass filter (BLPF). The ideal low pass filter has a very sharp transition like a circle while the gaussian and Butterworth filters have additional parameters that can make a softer transition as shown in the examples below.
def ddistance(u,v,p,q):
return (((u - (p/2)) ** 2) + ((v - (q/2)) ** 2)) ** (1/2)
# Ideal low pass filter of size (p,q)
def ideal_low_pass_filter(d0, p, q):
out = np.zeros((p,q))
for i in range(p):
for j in range(q):
out[i,j] = 1 if ddistance(i,j,p,q) <= d0 else 0
return out
# Gaussian low pass filter of size (p,q)
def gaussian_low_pass_filter(d0, p, q):
out = np.zeros((p,q))
for i in range(p):
for j in range(q):
out[i,j] = np.exp((-1*(ddistance(i,j,p,q) ** 2)) / (2*((d0) ** 2)))
return out
# Butterworth low pass filter of size (p,q)
def butterworth_low_pass_filter(d0, n, p, q):
out = np.zeros((p,q))
for i in range(p):
for j in range(q):
out[i,j] = (1) / (1 + (ddistance(i,j,p,q) / d0) ** (2*n))
return out
Examples of the 3 low pass filters: Ideal, Gaussian, and Butterworth with different parameters so you can see what the filters look like. You can see the the Ideal LPF has a very sharp change while the Gaussian and Butterworth LPF can be more gradual depending on what parameters are passed to the filters.
def generate_low_pass_filters():
ilpf_arr = [50, 100, 150]
glpf_arr = [10, 50, 100]
blpf_arr = [(50, 2), (50,4), (100,2)]
rows = 3
cols = 3
fig, axs = plt.subplots(rows, cols, figsize=(10,8))
for i, ax in enumerate(axs):
ax[0].imshow(ideal_low_pass_filter(ilpf_arr[i], 1000, 1000), cmap='gray')
ax[0].set_title(f'Ideal LPF: d0={ilpf_arr[i]}')
ax[1].imshow(gaussian_low_pass_filter(glpf_arr[i], 1000, 1000), cmap='gray')
ax[1].set_title(f'Gaussian LPF: d0={glpf_arr[i]}')
ax[2].imshow(butterworth_low_pass_filter(blpf_arr[i][0], blpf_arr[i][1], 1000, 1000), cmap='gray')
ax[2].set_title(f'Butterworth LPF: d0={blpf_arr[i][0]}, n={blpf_arr[i][1]}')
plt.suptitle(f'Low Pass Filters')
plt.tight_layout()
plt.savefig(f'low-pass-filters.png')
generate_low_pass_filters()
Down below I display the original images, the power spectrum and the log of the power spectrum of the FFT of the images. I also display the filter and apply that filter in the Fourier domain. I do the padding because that is what is described in the course textbook in section (4.7 The Basics of Filtering in the Frequency Domain, see Figure 4.35 and above it where they describe the process of filtering in the Fourier Domain if needed). It sounds like it is not mandatory but it can be used to apply different effects to the image. You can see that I organize the frequencies so that zero is in the middle of the FFT. You don't see much only looking at the power spectrum, but you do see more information when you look at the log of the power spectrum which has some strong lines. I also show the Gaussian low pass filter that is used to filter in Fourier domain and the output image after the filtering. You can see that since it removes some of the higher frequencies that the image becomes slightly blurred and more soft.
# Multiply image by (-1)^(x+y) to center the Fourier transform
def center_fourier_transform(img_in):
img = img_in.copy()
x, y = img.shape
for i in range(x):
for j in range(y):
img[i,j] = img[i,j] * ((-1) ** (i+j))
return img
def pad_image(img):
m,n = img.shape
padded_img = np.zeros((2*m,2*n))
padded_img[:m,:n] = img
return padded_img
def display_power_spectrum(img, axs):
# Display original image
ax = axs[0]
ax.imshow(img, cmap='gray')
ax.set_title('Original')
ax.set_xticks([])
ax.set_yticks([])
# Pad the image to size P and Q
m,n = img.shape
padded_img = pad_image(img)
# Center and computer fourier transform
fft = np.fft.fft2(center_fourier_transform(padded_img))
# Magnitude of the DFT is the frequency spectrum
# |F(u)| = [R(u)^2 + I(u)^2]^(1/2)
ps = ((fft.real ** 2) + (fft.imag ** 2)) ** (1/2)
# print(ps.shape)
# print(ps.max())
# print(ps.min())
# print()
# padded image
ax = axs[1]
ax.imshow(padded_img, cmap='gray')
ax.set_title('Padded Image')
ax.set_xticks([])
ax.set_yticks([])
# power spectrum, no log
ax = axs[2]
ax.imshow(ps, cmap='gray')
ax.set_title('Power Spectrum')
ax.set_xticks([])
ax.set_yticks([])
# power spectrum with log
# s = c*log(1+r)
ax = axs[3]
ax.imshow(np.log(1 + ps), cmap='gray')
ax.set_title('Log Power Spectrum')
ax.set_xticks([])
ax.set_yticks([])
# gaussian low pass filter
filter = gaussian_low_pass_filter(12, ps.shape[0], ps.shape[1])
ax = axs[4]
ax.imshow(filter, cmap='gray')
ax.set_title('Gaussian Low Pass Filter')
ax.set_xticks([])
ax.set_yticks([])
# Inverse FFT with filter
out = center_fourier_transform(np.fft.ifft2(fft * filter).real)
ax = axs[5]
ax.imshow(out, cmap='gray')
ax.set_title('Inverse FFT with Filter')
ax.set_xticks([])
ax.set_yticks([])
# Out Image
ax = axs[6]
ax.imshow(out[:m, :n], cmap='gray')
ax.set_title('Inverse FFT Image')
ax.set_xticks([])
ax.set_yticks([])
return True
def generate_power_spectrum(imgs):
rows = len(imgs)
cols = 7
fig, axs = plt.subplots(rows, cols, figsize=(16,14))
for i, ax in enumerate(axs):
idx = i
img = imgs[idx]
display_power_spectrum(img, ax)
plt.suptitle(f'FFT Power Spectrum')
plt.tight_layout()
plt.savefig(f'fft-power-spectrum.png')
generate_power_spectrum(images)
Down below you can see the process as described above except on the other images. You can see that some of the log power spectrums for some of these images have additional vertical and horizontal lines for some of the stronger lines in the image. There is also a few other smaller lines in the log power spectrum that somewhat mirror other lines in the spatial domain. The Gaussian low pass filter has the same impact of removing some of the higher frequencies and blurring the image once we convert it back to the spatial domain.
generate_power_spectrum(other_images[:5])
The below chart shows the impact of the three different low pass filters in the Fourier domain with different parameters on different images. Its interesting to note the slight ringing effect the top ideal LPF has on the image because of the sharp edge in the filter. The lower Gaussian filter really blurs the image while the higher parameters for the Gaussian filter just slightly blurs the image. It's the same for the Butterworth filter, depending on what parameters you pass to the function.
def apply_filter(img, filter):
# Pad the image to size P and Q
m,n = img.shape
padded_img = pad_image(img)
# Center and computer fourier transform
fft = np.fft.fft2(center_fourier_transform(padded_img))
# Magnitude of the DFT is the frequency spectrum
# |F(u)| = [R(u)^2 + I(u)^2]^(1/2)
ps = ((fft.real ** 2) + (fft.imag ** 2)) ** (1/2)
# Inverse FFT with filter
out = center_fourier_transform(np.fft.ifft2(fft * filter).real)
# Out Image
return out[:m, :n]
def generate_low_pass_filters_on_image(img):
ilpf_arr = [50, 100, 150]
glpf_arr = [10, 50, 100]
blpf_arr = [(50, 2), (50,4), (100,2)]
rows = 3
cols = 7
fig, axs = plt.subplots(rows, cols, figsize=(16,8))
for i, ax in enumerate(axs):
ax[0].imshow(img, cmap='gray')
ax[0].set_title('Original')
m,n = img.shape
p,q = m*2,n*2
filter = ideal_low_pass_filter(ilpf_arr[i], p, q)
ax[1].imshow(filter, cmap='gray')
ax[1].set_title(f'ILPF: d0={ilpf_arr[i]}')
ax[2].imshow(apply_filter(img, filter), cmap='gray')
ax[2].set_title(f'ILPF: d0={ilpf_arr[i]}')
filter = gaussian_low_pass_filter(glpf_arr[i], p, q)
ax[3].imshow(filter, cmap='gray')
ax[3].set_title(f'GLPF: d0={glpf_arr[i]}')
ax[4].imshow(apply_filter(img, filter), cmap='gray')
ax[4].set_title(f'GLPF: d0={glpf_arr[i]}')
filter = butterworth_low_pass_filter(blpf_arr[i][0], blpf_arr[i][1], p, q)
ax[5].imshow(filter, cmap='gray')
ax[5].set_title(f'BLPF: d0={blpf_arr[i][0]}, n={blpf_arr[i][1]}')
ax[6].imshow(apply_filter(img, filter), cmap='gray')
ax[6].set_title(f'BLPF: d0={blpf_arr[i][0]}, n={blpf_arr[i][1]}')
for j in range(cols):
ax[j].set_xticks([])
ax[j].set_yticks([])
plt.suptitle(f'Low Pass Filters On Image')
plt.tight_layout()
plt.savefig(f'low-pass-filters-on-image.png')
generate_low_pass_filters_on_image(images[0])
generate_low_pass_filters_on_image(other_images[0])
I implemented the phase correlation with the low pass filter. I show examples of the two original images, the phase correlation output, and what it looks like with different low pass filters and parameters. You can see in the raw phase correlation that there is usually a pixel that is brighter than the rest that shows the offset between the images. You can see that this offset actually becomes easier to see with some of the Ideal, Gaussian, and Butterworth low pass filters. This is very obvious for the last three where the bright pixel becomes more obvious to see with the LPF. On the top two, you can see the brightest pixel on the top left and then on the second one the brightest pixel is on the bottom right.
# Get the index of the max value
def get_max_index(pc):
max_idx = None
for i in range(pc.shape[0]):
for j in range(pc.shape[1]):
if max_idx is None:
max_idx = (i,j)
max_val = pc[i,j]
elif pc[i, j] >= max_val:
max_idx = (i,j)
max_val = pc[i,j]
return max_idx, max_val
def phase_correlation(f, g, ftype=0, filter_params={}):
m,n = f.shape
# Center for low pass filter
f_fft = np.fft.fft2(center_fourier_transform(f))
g_fft = np.fft.fft2(center_fourier_transform(g))
# G* = R(u) - jI(u)
g_star_fft = g_fft.real - g_fft.imag
a = f_fft * g_star_fft
# |F(u)|=[(R(u) ** 2) + (I(u) ** 2)] ** (1/2)
mag = np.sqrt((a.real ** 2) + (a.imag ** 2))
b = a/mag
# Use low pass filter
if ftype == 1:
b = b * ideal_low_pass_filter(filter_params['d0'], m, n)
elif ftype == 2:
b = b * gaussian_low_pass_filter(filter_params['d0'], m, n)
elif ftype == 3:
b = b * butterworth_low_pass_filter(filter_params['d0'], filter_params['n'], m, n)
return center_fourier_transform(np.fft.ifft2(b).real)
def generate_phase_correlation(img_pairs):
rows = len(img_pairs)
cols = 6
fig, axs = plt.subplots(rows, cols, figsize=(16,14))
for i, ax in enumerate(axs):
idx = i
f = img_pairs[idx][0]
g = img_pairs[idx][1]
ax[0].imshow(f, cmap='gray')
ax[0].set_title('f image')
ax[1].imshow(g, cmap='gray')
ax[1].set_title('g image')
ax[2].imshow(phase_correlation(f,g), cmap='gray')
ax[2].set_title('PC - no filter')
ax[3].imshow(phase_correlation(f,g,2,{'d0':50}), cmap='gray')
ax[3].set_title(f'PC - ILPF - d0=50')
ax[4].imshow(phase_correlation(f,g,2,{'d0':50}), cmap='gray')
ax[4].set_title(f'PC - GLPF - d0=50')
ax[5].imshow(phase_correlation(f,g,3,{'d0':50,'n':3}), cmap='gray')
ax[5].set_title(f'PC - BLPF - d0=50,n=3')
for j in range(cols):
ax[j].set_xticks([])
ax[j].set_yticks([])
plt.suptitle(f'Phase Correlation')
plt.tight_layout()
plt.savefig(f'phase-correlation.png')
img_pairs = [(images[4], images[5]), (images[1], images[0]), (other_images[0], other_images[1]), (other_images[2], other_images[3]), (other_images[0], other_images[3])]
generate_phase_correlation(img_pairs)
Down below I display a comparison on changing the parameters to the filters when doing the phase correlation and the output. You can see that as you change the parameters to the filters is really does highlight where the that peak in the phase correlation is, so there is probably an optimal low pass params for certain types of images when doing the filtering with the phase correlation.
def param_to_string(param):
# {'d0': 10, 'n': 2}
s = f'd0={param["d0"]}'
if param.get('n', '') != '':
s += f',n={param["n"]}'
return s
def generate_phase_correlation_params(f,g):
fnames = ['ILPF', 'GLPF', 'BLPF', 'BLPF_1', 'BLPF_2']
ftypes = [1, 2, 3, 3, 3]
fparams = {
'ILPF': [{'d0': 25},{'d0': 100},{'d0': 200}],
'GLPF': [{'d0': 25},{'d0': 50},{'d0': 100}],
'BLPF': [{'d0': 25, 'n': 1},{'d0': 50, 'n': 1},{'d0': 100, 'n': 1}],
'BLPF_1': [{'d0': 25, 'n': 3},{'d0': 50, 'n': 3},{'d0': 100, 'n': 3}],
'BLPF_2': [{'d0': 25, 'n': 5},{'d0': 50, 'n': 5},{'d0': 100, 'n': 5}],
}
rows = len(fnames)
cols = 6
fig, axs = plt.subplots(rows, cols, figsize=(16,14))
for i, ax in enumerate(axs):
idx = i
ax[0].imshow(f, cmap='gray')
ax[0].set_title('f image')
ax[1].imshow(g, cmap='gray')
ax[1].set_title('g image')
fname = fnames[i]
ftype = ftypes[i]
ax[2].imshow(phase_correlation(f,g), cmap='gray')
ax[2].set_title('PC - no filter')
fparam = fparams[fname][0]
ax[3].imshow(phase_correlation(f,g,ftype,fparam), cmap='gray')
ax[3].set_title(f'PC - {fname} - {param_to_string(fparam)}')
fparam = fparams[fname][1]
ax[4].imshow(phase_correlation(f,g,ftype,fparam), cmap='gray')
ax[4].set_title(f'PC - {fname} - {param_to_string(fparam)}')
fparam = fparams[fname][2]
ax[5].imshow(phase_correlation(f,g,ftype,fparam), cmap='gray')
ax[5].set_title(f'PC - {fname} - {param_to_string(fparam)}')
for j in range(cols):
ax[j].set_xticks([])
ax[j].set_yticks([])
plt.suptitle(f'Phase Correlation - Filter Params')
plt.tight_layout()
plt.savefig(f'phase-correlation-filter-params.png')
generate_phase_correlation_params(other_images[0], other_images[1])
I wrote out the functionality to perform the phase correlation and then find the pixel location with the maximal value. To then illustrate whether that implementation was correct I also align the images by resolving the four cases of ambiuity (as described in part 4 Mosaic Building). That way I can align the image using the peak found in the phase correlation and display the images together. You can see from the examples below that for my images using the pixel location with the maximal value works extremely well for aligning these images and it is very satisfying to see the images properly align. In my testing for my images, it seemed like using a simple threshold of 0.01 worked best for finding the peak value of the phase correlation and determining whether or not the images overlap. In practice, looking at the connected components from the phase correlation image thresholded at some value and then averaging the magnitude and position within those components in order to find the best peak to use for the two images works the best if you have correct thresholds. This is especially true depending on what low pass filter you use during the phase correlation and what parameters you pass to that filter. For example, if you look at the examples above you can see that certain filters make it so that the bright area in the phase correlation is a lot larger and consists of multiple pixels. In this scenario the base thing to do would be to use the thresholded connected components and do averaging on the position within those components find the best peak within that patch. Overall for most scenarios it seems like just using the pixel location with the maximal value performs very well as a simpler option to tune and use for these images.
import skimage
def pad_image_with_average(f,g):
f_temp = f.copy()
g_temp = g.copy()
fm, fn = f_temp.shape
gm, gn = g_temp.shape
if fm > gm:
temp = np.ones((fm,gn)) * g.mean()
temp[:gm, :gn] = g_temp
g_temp = temp
fm, fn = f_temp.shape
gm, gn = g_temp.shape
if gm > fm:
temp = np.ones((gm,fn)) * f.mean()
temp[:fm, :fn] = f_temp
f_temp = temp
fm, fn = f_temp.shape
gm, gn = g_temp.shape
if fn > gn:
temp = np.ones((gm,fn)) * g.mean()
temp[:gm, :gn] = g_temp
g_temp = temp
fm, fn = f_temp.shape
gm, gn = g_temp.shape
if gn > fn:
temp = np.ones((fm,gn)) * f.mean()
temp[:fm, :fn] = f_temp
f_temp = temp
return f_temp, g_temp
def normalized_correlation(a,b):
a = a - a.mean()
b = b - b.mean()
a_m = np.sqrt(np.square(a).sum())
b_m = np.sqrt(np.square(b).sum())
return ((a*b)/(a_m*b_m)).sum()
def find_best_corner(f, g, max_idx, display=True):
m, n = f.shape
gm, gn = g.shape
best = None
# Setup second image
out2 = np.zeros((m*2,n*2))
out2[max_idx[0]:max_idx[0]+gm, max_idx[1]:max_idx[1]+gn] = g
# Try the best
out1 = np.zeros((m*2,n*2))
out1[:m,:n] = f
out1[m:2*m,:n] = f
out1[:m,n:2*n] = f
out1[m:2*m,n:2*n] = f
tl1 = out1[max_idx[0]:m, max_idx[1]:n]
tr1 = out1[m:max_idx[0]+m, max_idx[1]:n]
bl1 = out1[max_idx[0]:m, n:max_idx[1]+n]
br1 = out1[m:max_idx[0]+m, n:max_idx[1]+n]
tl2 = out2[max_idx[0]:m, max_idx[1]:n]
tr2 = out2[m:max_idx[0]+m, max_idx[1]:n]
bl2 = out2[max_idx[0]:m, n:max_idx[1]+n]
br2 = out2[m:max_idx[0]+m, n:max_idx[1]+n]
tl = normalized_correlation(tl1,tl2)
tr = normalized_correlation(tr1,tr2)
bl = normalized_correlation(bl1,bl2)
br = normalized_correlation(br1,br2)
if display:
print('\n4 Patches Correlations:')
print(tl)
print(tr)
print(bl)
print(br)
arr = [tl, tr, bl, br]
best_idx = np.argmax(arr)
out = np.zeros((m*2,n*2))
if best_idx == 0:
out[:m,:n] = f
elif best_idx == 1:
out[m:2*m,:n] = f
elif best_idx == 2:
out[:m,n:2*n] = f
elif best_idx == 3:
out[m:2*m,n:2*n] = f
out[max_idx[0]:max_idx[0]+gm, max_idx[1]:max_idx[1]+gn] = g
best = out
return best, arr[best_idx], best_idx
def get_fitted_image(inp_img):
mask = (inp_img > 0).astype(int)
min_i = np.inf
min_j = np.inf
max_i = -1 * np.inf
max_j = -1 * np.inf
for i in range(inp_img.shape[0]):
for j in range(inp_img.shape[1]):
if mask[i,j] > 0:
if i < min_i:
min_i = i
if j < min_j:
min_j = j
if mask[i,j] > 0:
if i > max_i:
max_i = i
if j > max_j:
max_j = j
fitted_canvas = inp_img[min_i:max_i,min_j:max_j]
return fitted_canvas
def image_registration(f, g, display=True):
f_temp, g_temp = pad_image_with_average(f, g)
if display:
fig, axs = plt.subplots(1, 5, figsize=(14,4))
axs[0].imshow(f, cmap='gray')
axs[0].set_title('f image')
axs[1].imshow(g, cmap='gray')
axs[1].set_title('g image')
# Compute the phase correlation
pc = phase_correlation(f_temp, g_temp)
# Simple peak finding - max value index
max_idx, max_val = get_max_index(pc)
if display:
print(f'max_val:{max_val}')
if max_val < 0.01:
print(f'max val is less than threshold of 0.01, images dont overlap')
# Display phase correlation and max peak
if display:
# print(f'max: ({max_idx[0]},{max_idx[1]})')
axs[2].imshow(pc, cmap='gray')
axs[2].set_title('phase correlation with argmax')
# imshow flips an axis and plot doesnt
axs[2].plot(max_idx[1], max_idx[0], 'ro')
# Resolve the four cases of ambiguity after you find the peek
# based on the correlation on the overlapping patches in the
# spatial domain as described in Mosaic Building
best, best_corr, corner = find_best_corner(f, g, max_idx, display)
if display:
axs[3].imshow(best, cmap='gray')
axs[3].set_title('Aligned Images')
fitted_best = get_fitted_image(best)
if display:
axs[4].imshow(fitted_best, cmap='gray')
axs[4].set_title('Fitted Aligned Images')
if display:
for i in range(5):
axs[i].set_xticks([])
axs[i].set_yticks([])
plt.suptitle(f'Peak Finding - Argmax')
# plt.tight_layout()
plt.savefig(f'peak-finding-argmax.png')
return best, best_corr, max_idx, corner
def connected_components_peak_finding(pc, threshold=0.01):
print(pc.max())
pc_binary = (pc > threshold).astype(int)
print(f'num pixels greater than {threshold}: {pc_binary.sum()}')
# Connected components using the thresholded image
labels, num_labels = skimage.measure.label(pc_binary, 0, True)
print(f'num labels: {num_labels}')
# Get the region props
props = skimage.measure.regionprops(labels, pc)
print(f'len of props: {len(props)}')
# Round to get correct index
rounded = np.round(props[0].centroid_weighted)
# Using centroid weighted in order to get the centroid weighted with the intensity range.
return (int(rounded[0]),int(rounded[1]))
def image_registration_connected_components(f, g, display=True):
f_temp, g_temp = pad_image_with_average(f, g)
if display:
fig, axs = plt.subplots(1, 5, figsize=(14,4))
axs[0].imshow(f, cmap='gray')
axs[0].set_title('f image')
axs[1].imshow(g, cmap='gray')
axs[1].set_title('g image')
# Compute the phase correlation
pc = phase_correlation(f_temp, g_temp, 1, {'d0': 200})
# Connected components peak finding
# if not argmax_peak_finding:
max_idx = connected_components_peak_finding(pc)
max_val = pc[max_idx[0],max_idx[1]]
print(max_idx)
print(max_val)
if max_val < 0.01:
print(f'max val is less than threshold of 0.01, images dont overlap')
# Display phase correlation and max peak
if display:
# print(f'max: ({max_idx[0]},{max_idx[1]})')
axs[2].imshow(pc, cmap='gray')
axs[2].set_title('phase correlation with argmax')
# imshow flips an axis and plot doesnt
axs[2].plot(max_idx[1], max_idx[0], 'ro')
# Resolve the four cases of ambiguity after you find the peek
# based on the correlation on the overlapping patches in the
# spatial domain as described in Mosaic Building
best, best_corr, corner = find_best_corner(f, g, max_idx, display)
if display:
axs[3].imshow(best, cmap='gray')
axs[3].set_title('Aligned Images')
fitted_best = get_fitted_image(best)
if display:
axs[4].imshow(fitted_best, cmap='gray')
axs[4].set_title('Fitted Aligned Images')
if display:
for i in range(5):
axs[i].set_xticks([])
axs[i].set_yticks([])
plt.suptitle(f'Peak Finding - Connected Components')
# plt.tight_layout()
plt.savefig(f'peak-finding-connected-components.png')
return best, best_corr, max_idx, corner
Examples of peak finding by using the pixel location with the maximal value. I also printed and tested different thresholds for determining if the two images actually overlap. These examples below are specifically images that do overlap. I found that a value of 0.01 seemed to be a good threshold for highest pixel value for determining whether or not the images actually overlap. You can see the max values below which are between 0.02 - 0.26.
out = image_registration(images[4], images[5])
max_val:0.026673842529438815 4 Patches Correlations: 5.9705083025726226e-05 0.966170162052822 -0.1658247545475127 0.3542002832093809
out = image_registration(images[1], images[0])
max_val:0.05321337484149325 4 Patches Correlations: 0.0458111236335462 0.08873192701294352 0.962104436797337 0.07025478570448268
out = image_registration(other_images[1], other_images[4])
max_val:0.2252031832125647 4 Patches Correlations: 0.17668686030474332 0.03198615584644228 -0.05364424538843698 1.0
out = image_registration(images[1],images[4])
max_val:0.05893441811363937 4 Patches Correlations: 0.9684790749779867 0.3019129044452967 -0.05244891206039295 -0.06227593086926022
Examples of images that do not overlap for finding a good threshold of maximum value to be below a certain threshold. You can see that for these three examples that the max value is only up to around 0.007 which is under the threshold of 0.01. NOTE: I still display the images aligned improperly just to demonstrate that they do not overlap properly, but I print out that the images do not overlap based on that max val being below the threshold.
out = image_registration(images[1],other_images[3])
max_val:0.006266453972425221 max val is less than threshold of 0.01, images dont overlap 4 Patches Correlations: -0.006340864909351679 -0.029963400936686216 0.003736807406686117 -0.05463664066962649
out = image_registration(images[0],other_images[0])
max_val:0.0065026941749681145 max val is less than threshold of 0.01, images dont overlap 4 Patches Correlations: 0.03840708352365752 0.11247993084037902 -0.06313805591671455 0.0478074844440424
out = image_registration(images[4],other_images[2])
max_val:0.006687968061909436 max val is less than threshold of 0.01, images dont overlap 4 Patches Correlations: 0.11425482350366958 -0.0066857402146745515 -0.08631308961056362 0.054744001984065496
Examples of peak finding using connected components above a threshold of 0.01 and then averaging the magnitude and position within the found components in order to find the best peak for the offset between the two images.
out = image_registration_connected_components(images[4], images[5])
0.017702452726660628 num pixels greater than 0.01: 3 num labels: 1 len of props: 1 (107, 3) 0.017702452726660628 4 Patches Correlations: 5.9705083025726226e-05 0.966170162052822 -0.1658247545475127 0.3542002832093809
out = image_registration_connected_components(images[1],images[4])
0.04114481544165318 num pixels greater than 0.01: 5 num labels: 1 len of props: 1 (5, 360) 0.04114481544165318 4 Patches Correlations: 0.9684790749779867 0.3019129044452967 -0.05244891206039295 -0.06227593086926022
I created a function called mosaic_image_from_file() which will read in a list of images from a file, and then will create a mosaic from them and output them to a file. I also display the mosaic'd image in the notebook. The mosaic_image_from_file uses the mosaic_images function which computes the matrix of phase correlation and cross correlation from the aligned images and then uses that in order to align and put all images onto the canvas. You can see the mosaic created from the cell images and the beach trees images down below. In order to align the images properly, I use the technique in the last paragraph of (4 - Mosaic Building) on the assignment and resolve the four cases of ambiguity after the peak in phase correlation is found. I then computer the correlation on the overlapping patches in the spatial domain as described in the image on the homework. This helps resolve the four cases of ambiguity and figure out how to align the images.
def mosaic_images(imgs):
# compute matrix of cross correlation values
curr_images = imgs
out_imgs = []
corrs = np.zeros((len(curr_images), len(curr_images)))
offset_arr = []
corners = np.zeros((len(curr_images), len(curr_images)))
for i in range(len(curr_images)):
for j in range(len(curr_images)):
out, corr, offset, corner = image_registration(curr_images[i], curr_images[j], False)
out_imgs.append(out)
corrs[i,j] = corr
offset_arr.append(offset)
corners[i,j] = corner
print('matrix of cross correlation values for image pairs')
print(corrs)
# choose anchor image
anchor_idx = np.argmax(corrs.sum(axis=1))
print(f'anchor image idx: {anchor_idx}')
# put images on canvas using anchor image as base
m,n = curr_images[anchor_idx].shape
canvas = np.zeros((m*3,n*3))
so = m
canvas[so:so+m,so:so+n] = curr_images[anchor_idx]
for i in range(len(curr_images)):
if i != anchor_idx:
corner = corners[anchor_idx,i]
offset = offset_arr[(corners.shape[1]*anchor_idx)+i]
if corner == 0:
canvas[so+offset[0]:so+m+offset[0],so+offset[1]:so+n+offset[1]] = curr_images[i]
elif corner == 1:
canvas[so+offset[0]-m:so+m+offset[0]-m,so+offset[1]:so+n+offset[1]] = curr_images[i]
elif corner == 2:
canvas[so+offset[0]:so+m+offset[0],so+offset[1]-n:so+n+offset[1]-n] = curr_images[i]
elif corner == 3:
canvas[so+offset[0]-m:so+m+offset[0]-m,so+offset[1]-n:so+n+offset[1]-n] = curr_images[i]
fitted_canvas = get_fitted_image(canvas)
fig, axs = plt.subplots(1,2, figsize=(8,6))
axs[0].imshow(canvas, cmap='gray')
axs[0].set_title('Mosaic')
axs[1].imshow(fitted_canvas, cmap='gray')
axs[1].set_title('Fitted Mosaic')
plt.suptitle(f'Mosaic Building')
# plt.tight_layout()
plt.savefig(f'mosaic-building.png')
return canvas, fitted_canvas
def mosaic_image_from_file(path, fname, out_file):
imgs = []
with open(fname, 'r') as f:
for l in f:
ipath = path + l.strip()
imgs.append(convert_grayscale(plt.imread(ipath)))
mosaic, fitted_mosaic = mosaic_images(imgs)
plt.imsave(path + out_file, fitted_mosaic)
return
Below is an example of building the mosaic from images already read in an array in memory. You can see that it computes the cross correlation matrix and then uses that to select the anchor image and build the mosaic from there. I resolve the four cases of ambiguity as described above and in the assignment information.
out = mosaic_images(images)
matrix of cross correlation values for image pairs [[1. 0.96210444 0.06512203 0.97111869 0.97253864 0.05555162] [0.96210444 1. 0.95245673 0.97628675 0.96847907 0.96525974] [0.02836406 0.95245673 1. 0.10822695 0.98047127 0.92242006] [0.97111869 0.97628675 0.21037403 1. 0.94213967 0.0608977 ] [0.97253864 0.96847907 0.98047127 0.94213967 1. 0.96617016] [0.1322312 0.96525974 0.92242006 0.08987814 0.96617016 1. ]] anchor image idx: 4
out = mosaic_images(other_images)
C:\Users\arosX\AppData\Local\Temp\ipykernel_31980\992288798.py:39: RuntimeWarning: Mean of empty slice. a = a - a.mean() C:\Users\arosX\AppData\Local\Temp\ipykernel_31980\992288798.py:40: RuntimeWarning: Mean of empty slice. b = b - b.mean()
matrix of cross correlation values for image pairs [[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]] anchor image idx: 0
Below is the example of building a function that reads in a set of images specified in a file and builds and then saves a mosaic out to the desired location, based on the parameters that were passed into the function.
mosaic_image_from_file('./cell_images/', './cell_images/read_cells.txt', 'cell_mosaic.png')
C:\Users\arosX\AppData\Local\Temp\ipykernel_31980\992288798.py:39: RuntimeWarning: Mean of empty slice. a = a - a.mean() C:\Users\arosX\AppData\Local\Temp\ipykernel_31980\992288798.py:40: RuntimeWarning: Mean of empty slice. b = b - b.mean()
matrix of cross correlation values for image pairs [[1. 0.96210444 0.06512203 0.97111869 0.97253864 0.05555162] [0.96210444 1. 0.95245673 0.97628675 0.96847907 0.96525974] [0.02836406 0.95245673 1. 0.10822695 0.98047127 0.92242006] [0.97111869 0.97628675 0.21037403 1. 0.94213967 0.0608977 ] [0.97253864 0.96847907 0.98047127 0.94213967 1. 0.96617016] [0.1322312 0.96525974 0.92242006 0.08987814 0.96617016 1. ]] anchor image idx: 4
Below is the example of building a function that reads in a set of images specified in a file for the drone image and builds and then saves a mosaic out to the desired location, based on the parameters that were passed into the function.
mosaic_image_from_file('./other_images/', './other_images/drone.txt', 'drone_mosaic.png')
C:\Users\arosX\AppData\Local\Temp\ipykernel_31980\992288798.py:39: RuntimeWarning: Mean of empty slice. a = a - a.mean() C:\Users\arosX\AppData\Local\Temp\ipykernel_31980\992288798.py:40: RuntimeWarning: Mean of empty slice. b = b - b.mean()
matrix of cross correlation values for image pairs [[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.] [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]] anchor image idx: 0